- Кубический сплайн
-
Некоторая функция f(x) задана на отрезке , разбитом на части , . Кубическим сплайном дефекта 1 называется функция , которая:
- на каждом отрезке является многочленом степени не выше третьей;
- имеет непрерывные первую и вторую производные на всём отрезке ;
- в точках выполняется равенство , т. е. сплайн интерполирует функцию f в точках .
Для однозначного задания сплайна перечисленных условий недостаточно, для построения сплайна необходимо наложить какие-то дополнительные требования.
Естественным кубическим сплайном называется кубический сплайн, удовлетворяющий также граничным условиям вида:
Теорема: Для любой функции и любого разбиения отрезка cуществует ровно один естественный сплайн S(x), удовлетворяющий перечисленным выше условиям.Эта теорема является следствием более общей теоремы Шёнберга-Уитни об условиях существования интерполяционного сплайна.
Построение
Обозначим:
На каждом отрезке функция есть полином третьей степени , коэффициенты которого надо определить. Запишем для удобства в виде:
тогда
Условия непрерывности всех производных до второго порядка включительно записываются в виде
а условия интерполяции в видеОтсюда получаем формулы для вычисления коэффициентов сплайна:
Если учесть, что , то вычисление можно провести с помощью метода прогонки для трёхдиагональной матрицы.
Реализация на языке C++
#include <cstdlib> #include <cmath> #include <limits> class cubic_spline { private: // Структура, описывающая сплайн на каждом сегменте сетки struct spline_tuple { double a, b, c, d, x; }; spline_tuple *splines; // Сплайн std::size_t n; // Количество узлов сетки void free_mem(); // Освобождение памяти public: cubic_spline(); //конструктор ~cubic_spline(); //деструктор // Построение сплайна // x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены // y - значения функции в узлах сетки // n - количество узлов сетки void build_spline(const double *x, const double *y, std::size_t n); // Вычисление значения интерполированной функции в произвольной точке double f(double x) const; }; cubic_spline::cubic_spline() : splines(NULL) { } cubic_spline::~cubic_spline() { free_mem(); } void cubic_spline::build_spline(const double *x, const double *y, std::size_t n) { free_mem(); this->n = n; // Инициализация массива сплайнов splines = new spline_tuple[n]; for (std::size_t i = 0; i < n; ++i) { splines[i].x = x[i]; splines[i].a = y[i]; } splines[0].c = 0.; // Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки для трехдиагональных матриц // Вычисление прогоночных коэффициентов - прямой ход метода прогонки double *alpha = new double[n - 1]; double *beta = new double[n - 1]; double A, B, C, F, h_i, h_i1, z; alpha[0] = beta[0] = 0.; for (std::size_t i = 1; i < n - 1; ++i) { h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i]; A = h_i; C = 2. * (h_i + h_i1); B = h_i1; F = 6. * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) / h_i); z = (A * alpha[i - 1] + C); alpha[i] = -B / z; beta[i] = (F - A * beta[i - 1]) / z; } splines[n - 1].c = (F - A * beta[n - 2]) / (C + A * alpha[n - 2]); // Нахождение решения - обратный ход метода прогонки for (std::size_t i = n - 2; i > 0; --i) splines[i].c = alpha[i] * splines[i + 1].c + beta[i]; // Освобождение памяти, занимаемой прогоночными коэффициентами delete[] beta; delete[] alpha; // По известным коэффициентам c[i] находим значения b[i] и d[i] for (std::size_t i = n - 1; i > 0; --i) { double h_i = x[i] - x[i - 1]; splines[i].d = (splines[i].c - splines[i - 1].c) / h_i; splines[i].b = h_i * (2. * splines[i].c + splines[i - 1].c) / 6. + (y[i] - y[i - 1]) / h_i; } } double cubic_spline::f(double x) const { if (!splines) return std::numeric_limits<double>::quiet_NaN(); // Если сплайны ещё не построены - возвращаем NaN spline_tuple *s; if (x <= splines[0].x) // Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива s = splines + 1; else if (x >= splines[n - 1].x) // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива s = splines + n - 1; else // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива { std::size_t i = 0, j = n - 1; while (i + 1 < j) { std::size_t k = i + (j - i) / 2; if (x <= splines[k].x) j = k; else i = k; } s = splines + j; } double dx = (x - s->x); return s->a + (s->b + (s->c / 2. + s->d * dx / 6.) * dx) * dx; // Вычисляем значение сплайна в заданной точке по схеме Горнера (в принципе, "умный" компилятор применил бы схему Горнера сам, но ведь не все так умны, как кажутся) } void cubic_spline::free_mem() { delete[] splines; splines = NULL; }
Реализация на языке C# Платформа .NET// Интерполирование функций естественными кубическими сплайнами using System; class CubicSpline { SplineTuple[] splines; // Сплайн // Структура, описывающая сплайн на каждом сегменте сетки struct SplineTuple { public double a, b, c, d, x; } // Построение сплайна // x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены // y - значения функции в узлах сетки // n - количество узлов сетки public void BuildSpline(double[] x, double[] y, int n) { // Инициализация массива сплайнов splines = new SplineTuple[n]; for (int i = 0; i < n; ++i) { splines[i].x = x[i]; splines[i].a = y[i]; } splines[0].c = splines[n - 1].c = 0.0; // Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки для трехдиагональных матриц // Вычисление прогоночных коэффициентов - прямой ход метода прогонки double[] alpha = new double[n - 1]; double[] beta = new double[n - 1]; alpha[0] = beta[0] = 0.0; for (int i = 1; i < n - 1; ++i) { double h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i]; double A = h_i; double C = 2.0 * (h_i + h_i1); double B = h_i1; double F = 6.0 * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) / h_i); double z = (A * alpha[i - 1] + C); alpha[i] = -B / z; beta[i] = (F - A * beta[i - 1]) / z; } // Нахождение решения - обратный ход метода прогонки for (int i = n - 2; i > 0; --i) splines[i].c = alpha[i] * splines[i + 1].c + beta[i]; // Освобождение памяти, занимаемой прогоночными коэффициентами beta = null; alpha = null; // По известным коэффициентам c[i] находим значения b[i] и d[i] for (int i = n - 1; i > 0; --i) { double h_i = x[i] - x[i - 1]; splines[i].d = (splines[i].c - splines[i - 1].c) / h_i; splines[i].b = h_i * (2.0 * splines[i].c + splines[i - 1].c) / 6.0 + (y[i] - y[i - 1]) / h_i; } } // Вычисление значения интерполированной функции в произвольной точке public double Func(double x) { if (splines == null) return double.NaN; // Если сплайны ещё не построены - возвращаем NaN int n = splines.Length; SplineTuple s; if (x <= splines[0].x) // Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива s = splines[1]; else if (x >= splines[n - 1].x) // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива s = splines[n - 1]; else // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива { int i = 0, j = n - 1; while (i + 1 < j) { int k = i + (j - i) / 2; if (x <= splines[k].x) j = k; else i = k; } s = splines[j]; } double dx = (x - s.x); // Вычисляем значение сплайна в заданной точке по схеме Горнера (в принципе, "умный" компилятор применил бы схему Горнера сам, но ведь не все так умны, как кажутся) return s.a + (s.b + (s.c / 2.0 + s.d * dx / 6.0) * dx) * dx; } }
Литература
- Роджерс Д., Адамс Дж. Математические основы машинной графики. — М.: Мир, 2001. — ISBN 5-03-002143-4
- Костомаров Д. П., Фаворский А. П. Вводные лекции по численным методам.
- Волков Е. А. Глава 1. Приближение функций многочленами. § 11. Сплайны // Численные методы. — Учеб. пособие для вузов. — 2-е изд., испр.. — М.: Наука, 1987. — С. 63-68. — 248 с.
Ссылки
Категории:- Сплайны
- Численные методы
- Интерполяция
Wikimedia Foundation. 2010.